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We compare higher order gravity models to observational constraints from magnitude-redshift 
supernova data, distance to the last scattering surface of the CMB, and Baryon Acoustic Oscillations. 
We follow a recently proposed systematic approach to higher order gravity models based on minimal 
sets of curvature invariants, and select models that pass some physical acceptability conditions (free 
of ghost instabilities, real and positive propagation speeds, and free of separatrices). Models that 
satisfy these physical and observational constraints are found in this analysis and do provide fits to 
the data that are very close to those of the LCDM concordance model. However, we find that the 
limitation of the models considered here comes from the presence of superluminal mode propagations 
. for the constrained parameter space of the models. 
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CN ; I. INTRODUCTION 

For the last decade, several cosmological observations [l[ have been indicating that the expansion of the universe 



is accelerating. Cosmic acceleration constitutes one of the most challenging current problems in cosmology and all 
physics and several possible causes or explanations have been proposed, see for example the reviews [2|. These include 
i-C ■ various dark energy components, modifications to general relativity at cosmological scales, or a true or apparent effect 
of inhomogeneities in the universe. 

In this paper, we consider higher-order gravity models (HOG) that have attracted much attention in the recent 
literature since they can have a late- time self- acceleration as an alternative to a dark energy component [1, II], [H, B B HI • 
The models are derived from curvature invariants that are more general than the Ricci scalar and can exhibit early-time 
i inflation as well as late-time acceleration, e.g. 0, In these models, the dynamical equations are more elaborate 

\ than usual Fricdmann equations and the acceleration is a consequence of a different coupling between matter and 
t-H ■ space-time curvature. Some early papers, see e.g. [Til Il2^. have discussed the models within a more general context 
of unification theories for field quantization on curved back grou nds. From a point of view of phenomenology, some of 
the models were found to fit cosmological observations [l3l Il4, [H[ and other models fit solar system tests [l6| while 
other did not, e.g. [l7l |. 

ly-^ , However, in view of the large number of possible models, it became clear that a systematic approach to the models 
was very desirable 0, @]- Recently, the authors proposed in 0] a systematic approach based on minimal sets of 
invariants (bases) as the smallest independent buildin g bl ocks to construct these models. Indeed, previous theorems 
from the theory of invariants in general relativity [18. 1 19L |2G| showed that curvature invariants are not independent 
from each other and related by syzygies. For a given algebraic type of the Ricci tensor (e.g. the Segre classification 
[il], [23j]) and a given Petrov type of the Weyl tensor (i.e. symmetry classification of space-times) , e. o. [23l [24. [25Ll2q . 
there exists a complete minimal independent set (basis) of these invariants in terms of which all the other invariants 
may be expressed. Interestingly, all the Friedmann-Lemaitre- Robertson- Walker space-times of interest are covered by 
the Petrov Type O and Segre Type Al-[(1 1 1, 1)] (see for example [231 ]) where the number of independent invariants 
reduces to two invariants. 

In this paper, we use the systematic approach and select models free of ghosts, instabilities, and separatrices, and 
then compare the models to observations of supernova magnitude-redshift data, distance to the CMB surface, and 
Baryonic Acoustic Oscillations, including Hubble Key project and age constraints. 
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II. HIGHER-ORDER GRAVITY MODELS BASED ON MINIMAL SETS 

An approach for systematizing the study of HOG models was proposed by the authors, see [H|. The method is 
based on a connection made between theorems in invariants theory in relativity and higher-order cosmological models. 
We outline here the major ideas of this method and refer the reader to Q for the detail. The method involves the 
identification of the Petrov type corresponding to the symmetry of the spacetime, e.g. [13, 0, HE Eft] j and the Segre 
type of the Ricci tensor [2l|, [23]. In short, the Ricci tensor and the Weyl tensor are considered as tensorial operators 
acting on space-time vectors and bi-vectors, respectively, where the t ype s are determined from the multiplicity and 
properties of the eigenvalues and eigenvectors in each case, e.g. [23l. I24I l25l. l2t| . This limits the number of independent 
invariants that could be used to describe the spacetime. This unique number of invariants can be used to build a basis 
for describing the spacetimes in higher-order gravity models and thereby allowing for the development of a systematic 
(taxonomy) method of constraining the otherwise infinite number of invariants possible to consider. In fact, all 
Friedman-Lemaitre-Robertson- Walker (FLRW) manifolds fall under Petrov type O and Segre Type Al-[(1 1 1, 1)] 
with only two independent invariants, namely, {R, Rl} where 

R = g af3 R af 3, (1) 

is the Ricci scalar, and, the invariant 

m = ^s a f> s fl a , (2) 

is built using the mixed trace-free Ricci tensor, Sp a . 

Consistent with the basic idea, it was found in Q that these independent invariants eliminate some unwanted 
redundancies in the generalized Friedmann equations and lead to compact formulations. 

For illustration, let's give here the general field equations for f(R, Rl) models using the action of the form 

/ = M p [ d 4 x^g~[)-R + f(R,Rl)} + [ d 4 xV=g{L m +L r ), (3) 



*PJ - ~v ,L 2 - 

where M p = (— l/87rG), we vary the action with respect to the metric, g a p, and we find the generalized field equations 

S"? - \g aP R - \g aP f + f R S a P + \fn,g ap R + ff Q/3 /*; 7 1 ~ h, afi + \fnxS^S & 7 + \j RX S^R 
+\(f R iS^) n t + \g a(} (fm&% s - \{fRiS^\ a 7 - \{fm.Sn a ),^ = SttGT^, (4) 
where we have used the definitions 

f R= —f m = — (5) 

and T a @ is the energy-momentum tensor. 

In the field equations (HJ, the coupling between stress-energy and curvature is different than the one in GR (given 
by the first two terms on the LHS and energy momentum tensor on the RHS). The corresponding dynamics allows for 
some models to have late-time self-acceleration without the need for a cosmological constant or other forms of dark 
energy. In previous work Q, the authors were able to show a range of these models in the basis {R, Rl} that showed 
self-acceleration at late times using numerical and analytical methods. The discussion there focused on power-law 
type solutions with late-time acceleration free from separatrix singularities, so that the models can transition from a 
matter dominated phase to the desired cosmic acceleration phase seen by current observations. In the minimal set 
basis {R, Rl}, the authors of 4], built many models in a systematic way by following previous examples in which the 
spin-2 ghost instability was avoided by construction. We follow up here on that and expand the discussion to other 
instabilities. 



III. PHYSICAL CONSTRAINTS 



A ghost will refer to as a propagating degree of freedom which gives rise to a negative norm state upon quantization 
due to the wrong sign, leading to unphysical modes and particles [ill, HE US EH, [H, [33[ ■ As it was discussed in, 
for example [30| . while a necessary condition to avoid ghosts is that the decoupled equations be free of fourth order 
derivatives, the remaining second order derivatives must also have the correct sign in order to have a true ghost free 
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theory. In order to construct physically acceptable models, we restrict our models to the ones that comply with the 
two conditions and use some established results from previous studies [Tl|, [S3, HE HI, H3, H2, [H[. As usual, see for 
example [U H3, HI, Hil US, HI, Hf| , a theory with an action built from Ricci scalar plus a function of another invariant, 
for example the Gauss-Bonnet (GB) invariant, i.e. R + f(GB), can be re-written as the Einstein-Hilbert action plus 
the function coupled to a scalar field <j> with potential U(<fi), i.e. R + f(cf>) GB — U(4>). In the latter frame, the theory 
becomes like that of a Gauss-Bonnet one where the equations of motion (EOMs) then decouple into second order 
equations for the metric and for each scalar field involved. For example, (30j derived some conditions for higher order 
invariants so that the EOMs do decouple into second order equations for the metric and each scalar field. It is worth 
clarifying that the decoupled equations are not the background field equations in the R + f(GB) frame where the 
generalized Fricdmann equations contains higher order derivatives of the metric. For our present work, the conditions 
for ghost-free and non-superluminal propagations imposes some relations between the parameters that we are allowed 
to use in order to combine R and Rl. 

As discussed previously in Q, the fourth-order derivatives in the decoupled equations and the corresponding un- 
physical states are avoided by construction if we choose functions of the Gauss-Bonnet-like form 

ai R 2 + a 3 (--R 2 -8R1), (6) 
6 

where we have imposed the same condition 02 = — 4a3 used in [30l ] and where we used the vanishing of the Weyl 
tensor, see Next, to build models that avoid the separatrix, we choose to set a% ^ 03 which avoids the singularity 
at a — [4]. From now on, we limit the discussion to the case a\ ^ 123 = 1, and actions of the form 

m 4n+2 

^' m ) = -8*1]"' (7) 

that are free of fourth order derivatives by construction. 

Now, following the analysis of (3(j, we find, for our models, the equivalent equations to their equations (39), that 
are also necessary to have for a stable theory free of ghost, 

1 + |(6/3 - 8/ > 0, (8) 

l + ~(6/3-l)/# + 8ff/>0, (9) 

where H = a/a is the Hubble parameter and a dot stands for derivative with respect to cosmic time. 

Next, the spin-2 mode propagation speed is given by the ratio of the coefficients of the second derivatives, and the 
condition for a real positive and non-superluminal spin-2 propagation speed is thus given by 

, 1 + 1(6/3- l)fR + 8f , 

o<4 = — U < 1. (10) 

1 + 1(6/9- l)fR + 8Hf 

We also express, for our models, the condition for a real positive and non superluminal condition for the scalar 
mode propagation speed [3(| [U, [H, HH that reads 

r of i' i u\ 

(11) 



and 



2 32/# Kf-fH) 
< c 2 = 1 + —j- — j-i r- ^— - Yi~ hr- 



3(|(6/3 - l)(/i? + fR) + 8/iP) 3(1 + |(6/3 - l)fR + 8Hf) 

Now, we consider power-law solutions with a(t) oc t p in fiat FLRW spacetime for our accelerating branch as derived 
in our previous work Q and as given further by the solutions (I14p . For radiation dominated and matter dominated 
decelerating phases, we use p = 1/2 and p = 2/3, respectively. 

For the accelerating late time branch, we can write (30j . 

2 (6/8-l)p(2p-l) + 8(n+l)(4n + 3) 
C 2 Z7aa TTToI 1\ i oZTZ i i\ ' y Ll > 



and 



p(6/3-l)(2p-l)+8p(n+l) 
32p(n+l) 8(n+ l)(4n + 3-p) 



3[p(6/3- l)(2p- l)(4n + 2) + 8p 2 (n+ 1)] 3[p(6/3 - l)(2p - 1) + 8p(n + 1)] ' 



(13) 



4 



Inverse HOG 
Models 
for n > 


Self- Acceleration 
at late times. 
p > 1 


Deceleration During 
Matter Domination 
1/2 <p< 1 


Free 
of 

Separatrix 


Free of Ghost 
instability @ 


Free of Ghost 
instability @ 


Positive and real Spin-2 
propagation speed 
e| > 


Positive and real Sphi-0 
propagation speed 
cl >0 


P = 01 - § 
n = 1 


Oi < 0.621 


oi < 0.621 


Ol < | 


0.065 < oi < 0.621 


0.467 < oi < 0.621 


0.467 < oi < 0.621 


0.563 < oi < 0.621 


= oi - 1 
n = 2 


oi < 0.734 


oi < 0.734 


Ol < § 


0.474 < oi < 0.734 


0.636 < oi < 0.734 


0.636 < ai < 0.734 


0.683 < oi < 0.734 


/3 = 0l - | 
n = 10 


oi < 0.814 


Oi < 0.814 


Ol < | 


0.745 < oi < 0.814 


0.776 < oi < 0.814 


0.776 < oi < 0.814 


0.788 < oi < 0.814 



TABLE I: Physical parameter spaces for f(R,Rl) higher order gravity models (given by equation 0). Conditions for self- 
acceleration, deceleration during matter domination, a cosmic evolution free of separatrices, free of ghost instabilities, and 
requiring spin-2 and spin-0 modes to have propagation speed real and positive. The conditions are satisfied during the de- 
celeration and acceleration phases of the cosmic evolution (we use the late-time accelerating solution in equation (|14|l ; for 
radiation-dominated and matter-dominated decelerating phases, we use p = 1/2 and p = 2/3, respectively). For each con- 
straint indicated on the table, the allowed values for a\ are given. 



Finally, we recall here the discussion given in our previous work [4| for avoidance of separatrix singularities. The 
generalized Friedmann equations for actions of the form (J7J) were shown in to lead to the power law solutions with 

_ 1 y/6 12/3 + 3n + 42/3n + 48/3n 2 ± y/E 
P ~2 12^' 24/3(1 + n) ' ( } 



where 3 = 240/3n + 900/3V + 9n 2 + 588/3n 2 + 480/3n 3 + 2880/3V + 2304/3V + 24/3. (15) 

As discussed in in order to avoid singularities that lead to a separatrix, we must have no real solutions to the 
first two roots, (i.e. /3 < or cquivalcntly a\ < 5/6). The remaining solutions, a\ < 5/6, are divided into three 
regions depending on the type of attracting solutions. The first region gives negative attractors of no interest, and 
the second range has no real solutions, but the third range gives two real positive attractors with one large enough 
for self acceleration, see Q. 

Our results are summarized in Table I and include conditions for imposing a decelerating phase followed by an 
accelerating one, a ghost-free evolution, a real and positive propagation speed for all modes, and separatrix-free 
models. And we find models that meet all these conditions during both accelerating and decelerating phases. However, 
also when both phases are included in the analysis, we find that the parameter space range for which all the modes 
propagate subluminally is in conflict with either the separatrix-free condition or the ghost-free ranges. In other words, 
we find that if we impose on the models to satisfy all the conditions above during both deceleration and acceleration 
phases, then there are modes with superluminal propagations within one phase of the other of the cosmic evolution. 



IV. CONSTRAINTS FROM SUPERNOVA MAGNITUDE- REDSHIFT DATA, HUBBLE KEY 
PROJECT, AND BOUNDS ON AGE OF THE UNIVERSE 

One of the first compelling evidences for cosmic acceleration came from Supernovae type la observations. We use 
the recent Union data set of supernovae [48| to constrain the HOG models given in terms of the {R, Rl} basis. We 
use the distance modulus as a function of the redshift z, 

li[z) = m(z)-M = 51og 10 D L (z) + 25, (16) 

where m(z) is the magnitude of the supernova and M is considered as a nuissance parameter degenerate with the 
Hubble parameter, Hq, and also with our parameter rh defined below, Dl(z) is the luminosity distance in units of 
Mpc calculated for the HOG models and given by 

D L{Z ) = { 1 + Z )[^r )dz '. (17) 

where H(z') is the solution to the non-linear differential generalized Friedmann equation for HOG models (see for 
example Eq. (fig)) below) , similar to the Friedmann equation for GR models. We restrict ourselves to studying models 
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FIG. 1: LEFT: This figure shows how the approximation H approx of equation (|3U|) fits well the integrated H(z) from the ODE 
for our HOG models for redshifts > 5. RIGHT: The relative difference ~ g££ ^fe H shows that for z > 5 the fit is better than 
0.001. We use this approximation at high redshifts and a full numerical integration for z < 5. The relative difference between 
LCDM and HOG models is small at high redshift, and becomes important at late times. 



that fit the physical constraints discussed earlier in this paper, and in earlier work In view of the complexity of 
the analysis, we focus here on the cases of n = 1 and n = 2 in these models. 

For n = l, the generalized Friedmann equation reads (again, we use the notation (3 = ai — 5/6), 

3H 2 , ^- : ( 6048/3 2 H e H + 1152/5 2 iJ 8 - 2A0f3H 2 H 3 - 360f3H 4 H 2 

6(6/3# 2 + 2AI3H 2 H + 24/3ff 4 - H 2 f V 

-6H 2 H 3 + 5616f3 2 H' l H 2 + 3H 4 + 86A(3 2 H 5 H + 6HH 2 H + 1656(3 2 H 2 H 3 - !Uf3H 3 HH + 21Q0 2 HH 2 H 



-72{3HH 2 H + 864(3 H HH ~ 36/3i/ 4 + 108/3 2 £P + 48f3H°H + lAAflH Hj = 8nG Pm + 8nGp r . (18) 

We use in the analysis a modified version of the Markov Chain Monte Carlo (MCMC) package CosmoMC [4!|. We 
customized this package for our HOG models and use it to constrain their parameters. We found that an elaborate and 
essential step in our analysis was to numerically integrate the stiff ODE (|18[) in order to derive the Hubble expansion 
rate for a wide range of redshift. For this task, we followed the same notation and parameterizations as Ref. [l3l |. 
making a necessary change of variable in our integration and also using a good approximation at higher redshifts in 
order to obtain initial values for our codes to perform stable integrations from z = 5 down to z = 0. Indeed, in order 
to deal with the stiffness of the problem, we replaced H in (|T5]) by the logarithmic variable u = \n(H/rh) (see [l3j ) 
where 



Further, following |13| . we define 



and 



m = m/(12ai - 10) 1/6 . (19) 



12ai - 12 

= ~rz — , (20) 



12ai - 10' 



a = sgn(12 ai - 10), (21) 
which give the relationship for our parameter a\ as 

5a — 6 

ai = ^ TV 22 

o(a — 1) 
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FIG. 2: LEFT: HOG models with n = 1. 2D joint contour plots for (ai,Q m ) plane for, Union SNe la data sets, Hubble Key 
Project, and age of the universe where the inner and outer loops are 68% and 95%, respectively. CENTER: Same as LEFT, 
but for the (m,S7 m ) parameter space. RIGHT: ID parameter distributions for the same data set as at LEFT and CENTER. 



This allows one to write the sources term, see [la ], as 

S = ln (u r e- 4N +u m e~ 3N )/2, 

where N — In a and 

_ 8irG po 

3 171 

Further, with u> m — fl m h 2 and h = Hq/ (100km/ s/Mpc), one writes 

_ w m 
w m — To- 
rn 1 

Similarly, ui r is defined for radiation but we consider its contribution to be negligible at late times. 
Using this change of variable and the notation above, the stiff ODE (fT5|) can be re-written as [H 



u"yi(u')+y 2 (u') + 18a(y 3 (u')) 3 e 6u (e 



l\\3 6u( _2(u— u) 



1) 



0, 



where 



yi(x) = 8^27(oi - l) 2 x 2 + 18(6al - 5)(al - l)x + 2(9al - 7)(6al - 5)) j (6a a - 5) 2 
ya(i) = 4(l35(al - 1)V + 6(141al - 116) (al - l)x 3 + 2(6al - 5)(153al - 133)x 2 
+12(21al - 17)(6al - 5)x + 8(6al - 5) 2 ) j (6a x - 5) 2 
y 3 (i) = 2^3(0! - l)a; 2 + 2(6a x - 5)(x + 2)j j (6a x - 5), 



(23) 
(24) 

(25) 

(26) 

(27) 

(28) 
(29) 



and ' = d/dN. The final forms (|26)) - (|29|) above are similar to what was used in [l3| but expressed in terms of our 
parameters in the basis {Rl, R}. As mentioned earlier, unlike for equation (|18p . we found that the form above allows 
a stable numerical integration form redshifts of z — 5 down to z = 0. We also found it necessary to perform the 
integration forward in time (backward in the redshift) with initial conditions provided by the approximate solution 
given in [l3T ] and expressed in our notation as 



H, 



approx 



me" 1 + 



H" yi {u') + y 2 {u') 
36(7 (y 3 (u')) 3 



(30) 



As shown in Fig. 1, we verified that at higher redshifts (z > 2 — 3) the approximate solution provides an excellent fit 
to the numerical solution of the ODE (f!?6"]) . For our models, we thus use (|30|) in order to find the initial conditions for 
u and u! at z = 5 and then start a numerical integration of (|26p down to the supernova locations. Proceeding in this 
way, we obtained very stable programs to derive various Hubble plots for the HOG models. 
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FIG. 3: LEFT: HOG models with n = 2. 2D joint contour plots for (ai,Q m ) plane for, Union SNe la data sets, Hubble Key 
Project, and age of the universe where the inner and outer loops are 68% and 95%, respectively. CENTER: Same as LEFT, 
but for the (m,S7 m ) parameter space. RIGHT: ID parameter distributions for the same data set as at LEFT and CENTER. 



In a similar way, we derive the equations for the HOG models with n = 2. The modified Fricdmann equations read 

432 



u"y 4 (u')+y 5 (u') + 



(601 - 5) 



a(y 6 (u')fe 10u (e 2 ^ - 1) = 0, 



where 



y 4 (x) = 12(l5(ai - ifx 2 + 10(6al - 5)(al - l)x + 2(5al - 4)(6al - 5)) 
y 5 (» = 3(135(01 - 1) V + 36(23al - 19) (al - l)x 3 + 36(6al - 5)(8al - 7)x 2 
+8(27al - 22)(6al - 5)x + 4(6al - 5) 2 ) 
y 6 (x) = 3(oi - l)x 2 + 2(6oi - 5)(x + 1)) . 



And the expression for the approximation is given by 



H, 



approx 



3 fi (l + (6ai-5)- 



864a (y 6 {u r )) 4 



(31) 
(32) 

(33) 
(34) 

(35) 



We start here with Supernova data combined with the Hubble parameter value Hq = 72 ±8km/s/Mpc from Hubble 
Key Project [52j and the prior on the age of the universe so that only models with 10 < to < 20Gyr are considered. We 
use the Union set of supernovae which was compiled to gather the best supernovae from different surveys, including 
Supernovae Legacy Survey, ESSENCE Survey, HST, and other older sets [48|]. After selection cuts the 414 SNe la are 
reduced to 307. As usual, the fitting of the SNe la uses the standard \ 2 given by 



XSN 



i=307 

E 



(36) 



where a is the magnitude uncertainty and i the number of data points compared. 

Using the probes above, we use our modified version of CosmoMC to derive constraints on tt m , rh and ai for the 
HOG models with n = 1 and n = 2. 

For n = 1, we restrict the model parameter a\ to the region that passes the physical acceptability conditions given 
in Table I (see the last column). We find a best fit value for VL m = 0.396lg 242, m agreement with other studies using 
HOG models [10, E|, but in need of tighter constraints as we do in the next section. The best-fit value rh = 1.555 



+0.429 
-0.464 



is somewhat well constrained in the parameter space shown in Fig. 2, and within la deviations of twice the unitless 
Hubble constant h — Ho/lOOkm/s/Mpc, in consistency with [l3| for other models. While the fits do return a best 



value of ai = 0.575 



+0.046 
-0.012' 



we see in Fig. 2, the parameter a\ is not well constrained for this cosmological data. 



In fact, we find that all the physical parameter space for a\ fits well to the data. We find a best-fit Xsn ~ 311.7 
for our HOG model with fitting values given in Fig. 2. We procceed in a similar way for n = 2 models, and find 
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a best fit value O m = 0.239^o 62 (consistent with the one above within the error bars) and the model parameter 
rh = 1. 217^0 198 i s about ~ 2h again. The constraints are shown in Fig. 3. We find the best fit value for the model 
parameter a\ = 0.689lgQQg, but again, is not well constrained by the data presented in this section. All the physical 
parameter space allowed by the constraints in Table I fits well to the data with a Xsn ~ 311.4. The Xsn f° r n = 1 and 
n = 2 models are close to the Xsn ~ 308.1 obtained by the LCDM concordance model. In view of the of the possible 
systematic uncertainties in the supernova data, it is not clear that the difference between the two % 2 is significant. 
Our Hubble plot for the best fit HOG models with n = 1, 2 is presented in Fig. 4. 



V. ADDING CONSTRAINTS FROM DISTANCE TO THE CMB LAST SCATTERING SURFACE AND 

BARYON ACOUSTIC OSCILLATIONS 

In addition to supernova data constraints, Hubble Key project and age bounds, we consider here constraints from 
the distance to the CMB surface of last scattering and also Baryon Acoustic Oscillations (BAO). To do this, we use 
again a modified version of CosmoMC where various distance determinations are done using our numerical integrators 
as discussed above. We combine the approximation (|30|) at high redshift to the numerical integration of (F25|) at lower 
redshifts down to zero. 

Following [54|, we define three fitting parameters for comparison to WMAP5 data, the shift parameter, R, 

R(z*) = Vn^H Q (l + z*)D A (z*), (37) 
with the redshift, for the surface of last scattering, 

z* = 1048[1 + 0.00124(f2 b /i 2 )-°- 738 ][l + gi {Q m h 2 ) 92 ], (38) 

where 

0.0783(^ 2 )- - 238 
91 l + 39.5(f! fc /i 2 ) - 763 ' j 
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FIG. 5: LEFT: HOG models with n = 1. 2D joint contour plots for (clockwise from top left) (ai,fi m ), (m,£l m ), (rh,Ho), and 
(Hofim) planes for WMAP5, SDSS LRG(BAO), Union SNe la data sets, Hubble Key Project, and age of the universe where 
the inner and outer loops are 68% and 95%, respectively. On the top left, the region to the left of the vertical bar is the one 
satisfying the physical conditions of table I. RIGHT: ID parameter distributions for the same data set as at LEFT. The vertical 
line represents the cut for physical conditions in Table 1. 



and 

0.560 
1 + 21.1(SVi 

and third, the acoustic scale, l a , is 



•92 = ; I m 77* L2U . 8 i . (4°) 



l„ = (l + *)=^, (41) 

with the proper angular diameter distance, Da(z) — Dl(z)/(1 + z) 2 and the comoving sound horizon, r s (z„) 

r < z ) — / (42) 

VSJo a 2 H(a)^/T+WlUtth)a' 

with n 7 = 2.469 x lO" 5 /^ 2 for T cmb = 2.725K. 

Together the parameters Xi = (R,l a ,z*) are used to fit Xcmb = AxiCov -1 (xiXj)Axj with Axi — x, — x° bs and 
Cov~ x (xjXj)is the inverse covariance matrix for the parameters. 

Next, for the BAO, we follow [55[ and define the ratio of the sound horizon, r s (zd) to the effective distance, Dy as 
a fit for SDSS by 

2 _ ( r s (z d )/D v (z = 0.2) -0.198 x2 ( r s {z d )/D v {z = 0.35) - 0.1094 ^ 2 
Xbao ~ \ O0058 ) + V O0033 ) ' ^ 



with 



and the redshift, z& as 



where 



Dv{z) = {D\[z)(l + zf-^—\ , (44) 



-IsSp [1 + , ' ,W ' lfcl <45) 
h = 0.313(r> m /i 2 )~°- 419 [l + 0.607(r! m /i 2 ) - 674 ], (46) 
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FIG. 6: LEFT: HOG models with n = 2. 2D joint contour plots for (clockwise from top left) (ai,fi m ), (m,£l m ), (rh,Ho), and 
(Ho,Q m ) planes for WMAP5, SDSS LRG(BAO), Union SNe la data sets, Hubble Key Project, and age of the universe where 
the inner and outer loops are 68% and 95%, respectively. On the top left, the region to the left of the vertical bar is the one 
satisfying the physical conditions of table I. RIGHT: ID parameter distributions for the same data set as at LEFT. The vertical 
line represents the cut for physical conditions in Table 1. 



and 

b 2 = 0.238(a rn h 2 ) - 223 . (47) 

Our results for n=l are presented in Fig. 5 showing how the combination of SN+CMB_distance+BAO, along with 
Hubble Key Project and the age of the universe prior, provides tighter constraints on the HOG models considered. 
We find the best-fit Xsn+bao+cmb ~ 316.8 compared to Xsn+bao+cmb ~ 314.6 for the LCDM concordance 
model. We find Q m = 0.252±°; u2 ^ and to = 1.338±g^||. For the parameter a\, we obtain now better observational 
constraints with a\ — 0.5691q'oo6) but as shown on the Fig. 5a, the cuts there are still coming from the physical 
acceptability condition 0.563 < a\ < 0.621 (see Table 1). We also get tighter constraints on Hq — 71.00^ 1 '7Q. Next, 
for n = 2 models, we obtain a Xsn+bao+cmb ~ 316.2 and the model parameters are also better constrained with the 
additional data. The best fit values are f2 TO = 0.2531q;qis and to = 1.284t°;^. Also, the constraints on the model 
parameter a\ are now tighter, a\ = 0.692!°;^, but again the cuts come from the physical acceptability condition 
0.683 < a\ < 0.734. We also get tighter constraints on H = 70.551^21 ■ Our results for n = 2 arc in Fig. 6. Finally, 
we find again that because of the possible systematic uncertainties in the data considered, it is not clear that the 
difference between the \ 2 achieved by these models and the one from the LCDM model is significant. 



VI. CONCLUSION 

Following some initial work in we confront here some of the higher order gravity models (HOG), as expressed in 
terms of a minimal set of curvature invariants, to some physical acceptability conditions and observational constraints. 
The models studied had previously been constrained in [i[ by the requirement of a self-acceleration phase at late 
times, a deceleration phase during matter domination, and the avoidance of a separatrix region which would prevent 
a transition between the two evolutionary phases of the universe. Here, we add the further physical acceptability 
conditions for the models to be free of ghost instabilities and to have real and positive speeds of mode propagations. 
When these conditions are imposed, the allowable parameter space for the models considered is significantly reduced. 
We do find models that meet all these conditions during both accelerating and decelerating phases, but, the parameter 
space range for which all the modes propagate subluminally is found in conflict with either the separatrix-free condition 
or the ghost-free ranges. We find that if we impose on the models to satisfy all the conditions above during deceleration 
and acceleration phases, then they will have modes with superluminal propagations during one phase or the other. 
We compare these models to supernovae la, distance to the CMB surface, and BAO cosmological constraints. We 
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find a parameter space for the models considered that fit well the data and provide comparable fits to those achieved 
by the LCDM concordance model. The numerical integration over the full redshift range is elaborate for these 
models and required the use of methods to deal with stiff non-linear ODEs. We found that the combination of an 
approximate solution at high redshifts to a full numerical integration at low redshifts (z < 5) provides a stable and 
precise integration scheme to compare these and other HOG models to observations. While we find HOG models 
that pass the physical acceptability conditions above and cosmological constraints from supernovae, distance to CMB 
surface, and BAO constraints, it seems that the limitation of the models studied here comes from the presence of 
superluminal mode propagations. 
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